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Abstract 

Evolution of a background space-time metric and sub-horizon matter density perturbations in 
the Universe is numerically analyzed in viable f(R) models of present dark energy and cosmic 
acceleration. It is found that viable models generically exhibit recent crossing of the phantom 
boundary wde = — 1- Furthermore, it is shown that, as a consequence of the anomalous growth 
of density perturbations during the end of the matter-dominated stage, their growth index evolves 
non-monotonically with time and may even become negative temporarily. 
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I. INTRODUCTION 



The physical origin of the dark energy (DE) which is responsible for an accelerated ex- 
pansion of the current Universe is one of the largest mysteries not only in cosmology but 
also in fundamental physics [1]. Although the standard spatially flat A-Cold-Dark-Matter 
(ACDM) model is consistent with all kinds of current observational data [2], some tentative 
deviations from it have been reported recently [3, 4] which, if proven to be not due to system- 
atic and other errors, may eventually rule out an exact cosmological constant. Furthermore, 
in the ACDM model, the cosmological term is regarded as a new fundamental constant 
whose observed value is much smaller than any other energy scale known in physics. So, 
its understanding in fundamental physics is lacking today, although some non-perturbative 
effects may generate such a small quantity [5]. On the other hand, we know that "primor- 
dial DE," which is responsible for inflation in the early universe [6-8], is not identical to the 
cosmological constant, in particular, it is not stable and eternal. Hence it is natural to seek 
for non-stationary models of the current DE, too. 

Among them, f(R) gravity which modifies and generalizes the Einstein gravity by in- 
corporating a new phenomeno logical function of the Ricci scalar R, f(R), provides a self- 
consistent and non-trivial alternative to ACDM model, see e.g. Ref. [9] for a recent review. 
This theory is a special class of the scalar-tensor theory of gravity with the vanishing Brans- 
Dicke parameter oj bd [10, 11]. It contains a new scalar degree of freedom dubbed "scalaron" 
in Ref. [6] , thus, it is a non-perturbative generalization of the Einstein gravity. 

This additional degree of freedom imposes a number of conditions on viable functional 
forms of f(R). In particular, in order to have the correct Newtonian limit for R 3> Rq = 
R(to) ~ Hq where t is the present moment and H is the Hubble constant, as well as the 
standard matter-dominated stage with the scale factor behaviour a(t) oc t 2 ^ 3 driven by cold 
dark matter and baryons, the following conditions should be fulfilled: 

\f(R)-R\<&R, \f'(R) - 1| < 1, 1, J R>i? , (1) 

where the prime denotes the derivative with respect to the argument R. In addition, the 
stability condition f"(R) > has to be satisfied that guarantees that the standard matter- 
dominated Friedmann stage remains an attractor with respect to an open set of neighboring 
isotropic cosmological solutions in f(R) gravity. In quantum language, this condition means 
that scalaron is not a tachyon. Note that the other stability condition, f'(R) > 0, which 
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means that gravity is attractive and graviton is not a ghost, is automatically fulfilled in this 
regime. Specific functional forms that satisfy all these conditions have been proposed in 
Refs. [12-14] etc., and much work has been done on their cosmological consequences. 

In the previous paper [15] we calculated evolution of matter density fluctuations in viable 
f(R) models [12, 14] in the limiting case R ^> R during the matter-dominated stage and 
found an analytic expression for them. In this paper we extend the previous analysis and 
perform numerical calculations of the evolution of both background space-time and density 
fluctuations for the particular f(R) model of Ref. [14] without such restriction on R. As 
a result, we have found the phantom boundary crossing at an intermediate redshift z < 1 
for the background space-time metric and an anomalous behaviour of the growth index of 
fluctuations. 

The rest of the paper is organized as follows. In §2 we introduce evolution equations 
for the homogeneous and isotropic background and present results of numerical integra- 
tion. In §3 we report numerical solutions for the evolution of density fluctuations and other 
observables. Section 4 is devoted to conclusions and discussion. 

II. EVOLUTION OF THE BACKGROUND UNIVERSE 

We adopt the following action with a four-parameter family of f(R) models: 



where n, A, R s , and M are model parameters and S m is the action of the matter content 
which is assumed to be minimally coupled to gravity (thus, the action (2) is written in the 
Jordan frame). This is the model of Ref. [14] modified by the last term in (3) borrowed from 
the inflationary model of Ref. [6]. This term is introduced for several purposes associated 
with high-curvature behaviour of the theory. One of them, as explained in Ref. [14] , is 
to avoid excessive growth of the scalaron mass, m 2 s = 1/3 f"(R) in the regime (1), towards 
the early Universe, t — > 0. The other one is to remove the additional and undesirable 
"Big Boost" singularity which can arise in the original models [12-14] as was shown in Ref. 
[16] (see Refs. [15, 17, 18] for more discussion on this point). The value of M should be 




(2) 



(3) 
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sufficiently large in order not to destroy the standard cosmology of the present and early 
Universe. In particular, the values of M considered in Refs. [19, 20] are not high enough for 
this purpose, because M should not be smaller than the Hubble parameter H(t) during the 
iV ~ 60 last e-folds of inflation in the early Universe in order to avoid overproduction of relic 
scalarons, as well as to solve other cosmological problems. In fact, if we take M^3x 10 13 
GeV, the scalaron itself can act as an inflaton [6] and generate primordial scalar (adiabatic) 
and tensor perturbations [21, 22] with the amplitudes and slopes of their power spectra 
in agreement with all observational data available today. Note, however, that as shown in 
Ref. [18], such a "unified" model describing both primordial DE driving inflation in the early 
Universe and present DE driving recent acceleration of the Universe in the scope of f(R) 
gravity leads to slightly different predictions for parameters of the primordial perturbation 
spectra, as compared to the purely inflationary model with \R S = 0, due to a change in 
the number of observable e-folds of inflation N caused by different evolution of the Universe 
during generation and heating of usual matter after inflation. Furthermore, in this unified 
model the term in the square brackets in (3) should be modified for \R\ < R Q in such a way 
as to ensure the fulfillment of the stability condition f"(R) > in this region, too. 

So, we take this value of M and assume that the evolution of the Universe is identical to 
that in the standard ACDM model at high redshifts without any relic scalaron oscillations. 
Then the R 2 /QM 2 term is totally negligible in the epoch we are concerned here. Therefore, 
we do not include its contribution below. 

We can express field equations derived from the action in the following Einsteinian form. 

K - \KR = -SvrG (l? (m) + T? (DE) ) , (4) 

where 

8ttGT; (de) = F{R)R» - \HRK + (VV„ - 5»n)F(R), T(R) = f(R) - R (5) 

(the sign conventions here are the same as in Ref. [14]). Working in the spatially flat 
Friedmann- Robertson- Walker (FRW) space-time with the scale factor a(t), we find 

3H 2 = 8nGp - 3PH 2 + \{T'R - J 7 ) - 3Hp, (6) 

2H = -8nGp - 2T'H -P + Hp, (7) 

where H is the Hubble parameter and p is the energy density of the material content which 
we assume to consist of non-relativistic matter. 
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From (5) the effective energy density and pressure of dark energy can be expressed as 

87rGp DE = \{PR -F)- 3H 2 P - 3Hp' = -3HRP' + 3{H 2 + H)P - \l ', (8) 
87rG(p DE + P DE ) = 2HP - HP + P, (9) 

respectively, where R = 12H 2 + 6H. We define the DE equation of state parameter w DE by 
the ratio w DE = -Pde/pde- 

With the appropriate initial condition after cosmic inflation mentioned above, T takes 
an asymptotically constant value J 7 = —XR S at high redshift (apart from the R 2 /6M 2 
term which we neglect here). In this regime, evolution of the Universe is the same as that 
obtained from the Einstein action with a cosmological constant A(oo) = XR s /2. The scale 
factor therefore evolves as 



a = a '{^RT) smh, (v— y^-UJ • (10) 

where the suffix i denotes quantities at an initial time t = t{. 

The time dependence of p EE is mainly governed by the first term in the right-most 
expression of (8) initially. Since R < and T" > for stability, this means that the 
effective energy density of dark energy increases with time in this regime. Therefore, DE 
exhibits the phantom behaviour, w DE < — 1, during the matter-dominated stage with z > 1, 
which lasts only temporarily because the late-time asymptotic de Sitter stage has an effective 
cosmological constant smaller than A(oo). So, p EE stops growing after the end of the matter- 
dominated stage and begins to decrease. 

Indeed, as shown in Ref. [14], the late-time asymptotic de Sitter solution has a curvature 
R = Ri = x\R s where x\ is the maximal solution of the equation, 

*(l+* 2 ) w+1 (n) 
2[{l+x 2 Y+ 1 -l-{n + l)x 2 \ 1 J 

It satisfies the inequality x x < 2A, so that A(Rx) = i?i/4 < A(oo). These inequalities are 

saturated in the limit n ^> 1 for fixed x±, or x\ 1 for fixed n. In these cases cosmic 

evolution is indistinguishable from the standard ACDM model. 

Thus, this model naturally realizes crossing of the phantom boundary w DE = —1 in a 

recent epoch. Note that phantom behaviour of DE is generic in its models based on the 

scalar-tensor gravity [23] which includes the f(R) theory. Here we see that it is realized in 

all simplest stable f(R) models of present DE. 




n=2 A.=0.95 
n=3 A.=0.73 
n=4 X=0.61 
ACDM 




(a)Evolution of wde for different values of A with (b)Evolution of u>de for A m i n for n = 2,3, and 4. 
n = 2. 

FIG. 1: Evolution of the equation-of-state parameter of effective dark energy. 

The stability condition of this future de Sitter solution[24] , f'(Ri) > R\f"(R\), imposes 
the following constraint on x\. 



1 + x\) n+2 > 1 + (n + 2)x( + (n + l)(2n + l)x\, 



(12) 



which is stronger than any other constraint discussed above. For each n we can find x\ 
which marginally satisfies (12) and gives the minimal allowed value of A. Numerically we 
find (n, x lmin , A min ) = (2, 1.267, 0.9440), (3, 1.041, 0.7259), and (4, 0.9032, 0.6081) for each n, 
respectively (if n = 2, the analytic expression for Xi min is x\ mla = Vl3 — 2). For comparison, 
the analytic results for n — 1 are £i m i n = a/3 ~ 1.732, A min = 8/(3\/3) ~ 1.540. 

We numerically solve evolution equation (7) using (6) to check numerical accuracy, taking 
ti at the epoch when matter density parameter took Qi = 16TrGpi/(16irGpi + XR S ) = 0.998. 
We determine the current epoch by the requirement that the value of Q takes the ob- 
served central value Q$ = 0.27 and R s is fixed so that the current Hubble parameter 
H = 72km/s/Mpc is reproduced. We find the ratio R s /Hq is well fit by a simple power- 
law R s /Hl = c n \- pn with (n,c n ,p n ) = (2,4.16,0.953), (3,4.12,0.837), and (4,4.74,0.702), 
respectively, whereas in the ACDM limit it would behave as R s / Hq = 6(1 — Oq) /A ~ 4.38A -1 . 

Figures 1 depict evolution of u>de as a function of redshift z where phantom crossing is 
manifest. As expected, it approaches wbe = — 1 = constant as we increase A for fixed n. 
For minimal allowed values of A, deviations from wde = — 1 are observed at ~ 5% level in 
both directions for z < 2 independently of n. Such behaviour of »de is well admitted by all 
most recent observational data, see e.g. Ref. [2]. The average value of wde over the interval 



< z < 1 to which all BAO and most of SN data refer is very close to —1. Moreover, in 
this range (but not for larger values of z), the behaviour of wde for minimal allowed values 
of A (i.e. for largest possible deviations from the ACDM background model) is well fitted by 
the CPL fit[25] w DE {z) = w + w a z/{l + z) with (n, A min , w , w a ) = (2,0.95,-0.92,-0.23), 
(3, 0.73, -0.94, -0.22) and (4, 0.61, -0.96, -0.21), respectively. |1 + w \ and \w a \ decrease 
slowly for larger values of n. These values of wo and w a lie very close to the center of the 
68% and 95% CL ellipses for all combined data in Fig. 13 of Ref. [2]. 

As explained above, this phantom crossing behaviour is not peculiar to the specific choice 
of the function (3) but a generic one in models which satisfy the stability condition T" > 0. 
Indeed, a similar behaviour has been observed in other f(R) DE models, too [18, 26]. We 
also note that different definitions of pde, Pbe, and wde have been used in literature [27] 
which lead to different behaviour of u>de- 

Although the behaviour of dark energy is quite different depending on model parameters, 
the total expansion factor a /ai from the epoch f2j = 0.998 to the present varies only between 
ao/aj = 10.8 and 11, the latter corresponding to the value in the ACDM model. 

We have also calculated the quantity B(z) = (/" / f')(dR/d\nH) introduced in Ref. [28] 
at present time. We have found B{0) = 0.21, 6.1 x 10~ 5 , and 0.17, for (n, A) = (2,0.95), 
(2,8), and (4,0.61), respectively. 



III. DENSITY FLUCTUATIONS 



We now turn to evolution of density fluctuations. In f(R) gravity, the evolution equation 
of density fluctuations, 5, deeply in the sub-horizon regime is given by [29, 30] 

5 + 2H5 - AitG&pd = 0, (13) 

where 

n i + 4^EL 

G cS = p 1± ^ , F(R) = f(R). (14) 

This equation reduces to the correct evolution equation for all wavenumbers for the CDM 
model in the Einstein gravity where F — 1. 

In the previous paper [15] we obtained an analytic solution in the high-curvature regime 
when the scale factor evolves as a(t) oc t 2 / 3 and F takes the asymptotic form 

\ -2n-l / jj, \ -JV-1 



with the following correspondence: 



N = 2n and R c = R s (2n\) 



l/(2n+l) 



(16) 



The two independent solutions of (13) in this regime read 



S k (t) 



3 (ji±M ( lT + * ,3 \ ()7) 



/ ±5- y^33 ±5 + y/33 5 
X 2 1 ( v 4(3iV + 4) , 4(3iV + 4)' 2(3iV + 4)' " a 2 R 2 \t 

in terms of the hypergeometric function [15]. In the following discussion, we consider the 
upper sign solution only, because the other solution corresponds to the decaying mode and 
is singular at t — > 0. Then the solution behaves as 

-l+y/33 



2 

S k (t) ^ 8 ik ( I) 3 and 5 k (t) ^ 8 ik C(k) 



respectively. The transfer function, C(k), is given by 



C(k) 



r i + 



2(3iV+4) 



33 \ 



2(3^+4) I 



1 l^ 1 "'" 4(3V+4) y 1 ^4(3iV+4) 



3(iV+ l)fc 2 /3,R c t 2 



a 2 i? c 



2 \ ^+2' 



(18) 



-5+\/a3 

4(3AT+4) 



r i + 



4(3n+2) 



r 



33 



4(3n+2) 



( -i , _5+yj3_A r / _5W33_\ 
\ 8(3n+2) y 1 ^8(3n+2) y 



6nA(2n + 1)A; 2 /3i? s t 



am. 



2\ 2(n+l)' 



8(3n+2) 



where 



Note that the effective gravitational constant (14) reads 



GW = G 1 + - 



1 k 2 /a 2 m 2 s 



(19) 



(20) 



(21) 



3 1 + k 2 /a 2 m 2 / 

in the high-curvature regime when F = 1. In the position space, such a theory has the 
potential 

V(r) = "f (l + , (22) 

per unit mass [33] for such sufficiently small r for which time dependence of m s (t) may be 
neglected. Thus, each Fourier mode feels 4/3 times the conventional gravitational force if 
and only if k/a(t) > m s (t) = (3F')~ 1/2 . 
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The transition from former temporal behaviour to the latter one in (18) occurs at the 
epoch tk determined by 

n+l 



k = a(t k )m s {tk) = a(tk) 



R(t h 



(23) 



Qn(2n + l)\J \ R s 

The above expression is proportional to t k 2n 4 ^ 3 for those modes which physical wavenumber 
(momentum) k/a(t) crosses the scalaron mass m s (t) in the high-curvature regime. This 
explains /^-dependence of the transfer function (19) [14]. If we adopt an expression of R(t) 
in ACDM, 



R(t) 



3H 2 



n 



mO 



a(t) 



+ 4(l-ft m0 ) 



(24) 



we can further approximately obtain the crossing time, t*(k), for a smaller wavenumber, fc*, 
as well: 



y/Gn(2n + l)£ + * 



3Q 



mO 



a(t*) 



+ 12(1 - a rn0 ) 



n+l 



H . 



(25) 



From (25) we find that the physical wavenumber crossing the scalaron mass today is given 
by 

' 3.2\ im H n 



k 9.57 n+1 A( n+ ^ p ""^ 



fl o v/6n(2n + l)cZ 



n+ 



5.3X 2A3 H 
5.0\ 2m H 



(n = 2) 
(n = 3) . 
(n = 4) 



(26) 



Thus, except for cases with large A, all observable scale feels the scalaron force today. 

Since the analytic solution (17) is valid in the high-curvature era only, we must solve (13) 
numerically to obtain a full solution using the analytic solution as an initial condition. Figure 
2 depicts the ratio of linear density fluctuation in f(R) model, S^q, to that in the ACDM 
model, <5acdM) with the same initial condition. Fluctuations with small wavenumbers have 
practically the same value as those in the ACDM model, while those on larger wavenumbers 
acquire additional growth due to the scalaron force with the additional power k 4< - 3n+2 '> as 
given in (19). From (26), the physical wavenumber of this transition is given by 



a 



1.07 x 10- 3 /iMpc _1 
8.44 x 10- 3 /iMpc _1 
8.12 x lO^/iMpc -1 



(n = 2, A = 1) 
(n = 2, A = 3) 
(n = 2, A = 10), 



(27) 



that explains the figure well. 
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FIG. 2: The ratio of linear density perturbations 5fRc/^ACDM at present as a function of k for 
three different values of A with n = 2. 

In order to make a simple comparison of our results with observations of galaxy clustering, 
we define an effective wavenumber, k e s(r), corresponding to each length scale r, in terms of 
the top-hat mass fluctuation within the same radius: 



ol 



d " k -\W{kr)\ 2 P{k) = ^sfP (k eS (r)) , W(kr) = 3jl(kr) 



(28) 



(2tt) 3 ' v n r/ "(2f) 3 * v * v/ " kr 

Here P(k) is the linear matter spectrum obtained by the standard CDM transfer functional] 
with the scale-invariant initial power spectrum of perturbations, i.e. with the primordial 
spectral index n s = 1, and W(kr) is the Fourier transform of the top-hat window function. 

The wavenumber of our particular interest is the scale corresponding to <jg normalization, 
for which we find k e s(r = 8/i _1 Mpc) = 0.174/iMpc -1 . Figure 3 depicts the redshift evolution 
of the ratio ^fRc/^ACDM for this scale for the same values of n and A as in Fig. 2. Note that 
this ratio does not stop growing at the accelerated stage of the Universe expansion which 
begins at z ~ 0.8 for A = 1 and z ~ 0.75 for two other values of A. Since the standard ACDM 
model normalized by large-scale CMB observations explains galaxy clustering at small scales 
well, <5fRG should not be too much larger than <5acdm at these scales. We may typically require 
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FIG. 3: The ratio of linear density perturbations <5fRG / <5acdm(& = 0.174/iMpc ) as a function of 
redshift for three different values of A with n = 2. 

(^ircMacdm) 2 ^ = 0.174/iMpc^ 1 ) < 1.1. Although we neglect non-linear effects here, the 
difference between linear calculation and non-linear N-body simulation remained smaller 
than 5% at the wavenumber 0.174/iMpc~ 1 [32]. 

Figures 4 represent (^fRc/^ACDM) 2 ^ = 0.174/iMpc -1 ) as a function of A for n = 2, 3, and 
4. From the analytic formula (19), this A dependence would have the form (tffRc/^ACDM) 2 <x - 

(2n- Pn + l)(V33-5) 

C (k) oc A 4(3n+2) which is depicted by a broken line in each figure. This curve, 
however, does not match the asymptotic behaviour {Strq/ $acdm) 2 — > 1 for large A. We 
find that an exponential function 

( W^acdm) 2 = 1 + Ke^ x (29) 

fits the numerical calculation very well with (n,b n ,q n ) = (2,0.47,0.19), (3,0.43,0.49), and 
(4, 0.39, 0.70) , respectively. From these figures, in order to keep deviation from ACDM 
model smaller than 10% at k — 0.174/iMpc -1 , we find A should be larger than 8.2, 3.0, and 
1.9 for n = 2, 3, and 4, respectively. 

From these analysis, we can constrain the parameter space as Fig. 5. The region which 
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k = 0.174 hMpc"' 
numerical fit 
analytical fit 




k = 0.174 hMpc 
numerical fit 
analytical fit 




(a)n 



(b)n = 3. 



k = 0.174 hMpc 
numerical fit 
analytical fit 




(c)n = 4. 

FIG. 4: The present ratio (5fRG/^ACDM) 2 (^ = 0.174/iMpc -1 ) as a function of A together with two 
fitting functions. 



satisfy (8^/5. 



acdm) 2 (& = 0.174/iMpcT 1 ) < 1.1 corresponds to above the solid line. We also 
show the 20% boundary by the broken line. The region below the dotted line is forbidden 
because of instability of the de Sitter regime. 

Next we turn to another important quantity used to distinguish different theories of 
gravity, namely, the gravitational growth index, 7(2), of density fluctuations [33-38]. It is 
defined through 



din 5 



n m (zy (z \ or 7 (*) V 



din a " Ly ' ' lx ' logf2 m ^ ^ 

It takes a practically constant value 7 = 0.55 in the standard ACDM model[34], but it evolves 
in time in modified gravity theories in general. We also note that j(z) has a nontrivial k- 
dependence in f(R) gravity since density fluctuations with different wavenumbers evolve 
differently. Therefore, this quantity is a useful measure to distinguish modified gravity from 
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FIG. 5: Constraints for parameter space. 



the ACDM model in the Einstein gravity. 

Figures 6 show evolution of 7(2) together with that of G b q/G for different values of k. In 
the early high-redshift regime, j(z) takes a constant value identical to the ACDM model be- 
cause f(R) gravity is indistinguishable from the Einstein gravity plus a positive cosmo logical 
constant then. It gradually decreases in time, reaches a minimum, and then increase again 
towards the present epoch. We can understand this tendency from the evolution equation 
for 7(2) [36], 



[1 + z) In(l-n DB )^ 

= -(1 - Q DE y - ~ [1 + 3(2 7 - 1 Vde^de] + ~(1 



n 



DE 



,1-7 



(31) 



where Qde — 1 — is the density parameter of dark energy based on (8). In the high- 
redshift era when f2 DE is small, the above equation may be approximated as 

c?7 



[1 + z)n DE 

3 / G c s 
~2\~G~ 



dz 



1 +n 



DE 



11 

T 



7 



-V-d 

11 ) 2 y 



7) 



G 



cff 



G 



-(2 7 -l)( WDE +l) 



(32) 
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(a)Evolution of 7(2) for n = 2 A = 1. 
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(b)Evolution of G eff /G for n = 2 A = 1. 
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(c)Evolution of 7(2) for n = 2 A = 3. (d)Evolution of G cff /G for n = 2 A = 3. 

FIG. 6: Evolution of 7(2) and G cS /G. 

In the earlier stage, the first term in the right-hand side is more important. That explains 
why 7(2) starts to decrease when G e e/G starts to increase. As time goes by towards lower 
redshifts, the second term becomes more important to make 7(2) increase again. We note 
that recently Narikawa and Yamamoto[38] calculated time evolution of 7(2) in a simplified 
model (15) numerically and also obtained some analytic expansion, which behaves qualita- 
tively the same as our numerical results but with much more exaggerated amplitudes. Our 
results, which satisfy all the viability conditions, exhibit milder deviation from the ACDM 
model than those they found. Existing constraints on the growth index[39] are not strong 
enough to detect any deviation from the ACDM model and/or to obtain new bounds on 
f{R) DE models, but future observations may reveal its time and wavenumber dependence. 

Another quantity which can characterize the evolution of density perturbations more 
directly is the ratio 5[rg( z = 0-5)/8{rg( z — 0). However, it varies only from 0.75 to 0.78 for 
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FIG. 7: The evolution of = with n = 2. 

different choices of the model parameters when the current matter density parameter is fixed 
to fl m Q = 0.27 and n > 2. This variation is smaller than that caused by the uncertainty of 
^mo[33]. So, at present it does not help much to single out the best DE model among the 
considered ones, in contrast to the f(R) DE model[12] (it has the same behaviour (15) for 
R 3> R s ) in the case corresponding to n = 0.5 in our notations which was recently studied 
in Ref. [40]. 

Finally we consider the quantity I/77 = namely the ratio of gravitational potential 

to curvature perturbation, for which some results from observational data were recently 
obtained in Ref. [4]. In f(R) gravity, 1/rj is expressed as 



2 - 



V 



1 1 9 k 2 F' ' 

a 2 F 



(33) 



Due to the stability conditions F > 0, F' > 0, this quantity always lies between 1 and 2. 
Thus, stable f{R) DE models may not explain such a large value of 1/77 which is presented 
in Ref. [4] for the redshift interval 1 < z < 2. Figure 7 shows the evolution of 1/rj for n = 2 
and A = 0.95 (the minimal possible value) and 8. 
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IV. CONCLUSIONS 



In the present paper we have numerically calculated the evolution of both homogeneous 
background and density fluctuations in a viable f(R) DE model based on the specific func- 
tional form proposed in Ref. [14]. We have found that viable f(R) gravity models of present 
DE and accelerated expansion of the Universe generically exhibit phantom behaviour during 
the matter-dominated stage with crossing of the phantom boundary wde = — 1 at redshifts 
z < 1. The predicted time evolution of wbe has qualitatively the same behaviour as that 
was recently obtained from observational data in Ref. [3] . However, it is important that 
the condition of stability, or even metastability, of the future de Sitter epoch strongly re- 
stricts possible deviation of w DE from —1 by several percents in these models. Thus, the DE 
phantomness should be small, if exists at all, that agrees with present observational data. 
Still for the models considered, it is not so hopelessly small as in the case of the similar 
model[12] with n = 0.5 recently considered in Ref. [40] using data on cluster abundance. 
Note also, that in contrast to Ref. [41], we do not impose the so called thin-shell condition 
\A(f'(R) — 1)| < |$at|, where & N is the Newtonian potential of matter inhomogeneities 
and A means change in the quantity in question, for scales exceeding galatic ones where 
a background matter density approaches the cosmological one. On the other hand, this 
condition is satisfied automatically for matter overdensities more than 10 for the parameter 
range n > 2 considered in our paper. 

As for the density fluctuations, we have numerically confirmed our previous analytic 
results of a shift in the power spectrum index for larger wavenumbers which exceed the 
scalaron mass during the matter-dominated epoch[15], while for smaller wavenumbers fluc- 
tuations have the same amplitude as in the ACDM model. Once more, the future de Sitter 
epoch stability condition bounds possible increase in density fluctuations for cluster scales 
(compared to the ACDM model) by ~ 20% for n > 2. On the contrary, if it is proven from 
observational data that this increase is less than 5%, then the background evolution should 
be practically indistinguishable from the ACDM one: \w-de + 1| < 10~ 4 for n = 2. This 
shows that a 8 and related density perturbations tests are the most critical ones for the f(R) 
DE models considered in the paper. We have obtained that the upper limit on |u>de + 1| for 
n = 2 and A = 8 is 4.4 x 10~ 5 when z = 0.16, which is of the same order as -B(O). 

We have also investigated the growth index j(k, z) of density fluctuations and have pre- 



16 



sented an explanation of its anomalous evolution in terms of time dependence of G e g. Since 
7 has characteristic time and wavenumber dependence, future detailed observations may 
yield useful information on the validity of f(R) gravity through this quantity, although 
current constraints have been obtained assuming that it is constant both in time and in 
wavenumber [4, 39]. Another related observational test of this model is supplied by the 
large-scale structure of the Universe which should be different from that in the ACDM 
model. In particular, voids are expected to be more pronounced since the effective gravi- 
tational constant is bigger inside them compared to large matter overdensities where it is 
practically equal to that measured in laboratory. 
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